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temperature phase boundaries without requiring the calculation of the whole free 
energy surface of the alloy system. 
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1. Introduction 

Monte Carlo (MC) simulations [0-^ of lattice models are a widely used method to 
compute thermodynamic properties of substitutional alloys. A lattice model MC 
simulation only requires, as an input, a so-called cluster expansion which 
defines the energetics of the system by specifying the energy associated with any 
atomic arrangement on a given lattice. Cluster expansions are typically constructed 
from a fit to the energies of a set of structures calculated, for instance, from first- 
principles density-functional-theory calculations P], thus enabling the prediction of 
thermodynamic properties even in the absence of experimental data. 

In its simplest form, a cluster expansion specifies the energy difference between 
nearest-neighbor chemical bonds joining identical and distinct species in a binary alloy 
{Eab — {Eaa + Ebb) /2) and predicts compound formation energies from the number 
of each type of nearest-neighbor bond. However, this formalism can be extended to 
allow for arbitrary interaction ranges and multibody interactions so that the energetics 
of the alloy can be modeled with an arbitrary accuracy by including a sufficient number 
of such interactions. It has been demonstrated that the energetic contributions arising 
from displacements of the atoms away from their ideal lattice sites P,|7|,p|-p!^ can be 
accounted for within this framework. Moreover, by allowing temperature-dependent 
interactions, entropic contributions arising from vibrational and electronic excitations 
can be included as well (see for a review and |T^-|TP[ for recent examples). Since 



long-range interactions are often needed to properly describe the energetics of an alloy 
system, MC simulations have become a preferred method to compute thermodynamics 
properties from cluster expansions, because it would be difficult to account for such long- 
range interactions within accurate mean-field methods, such as the Cluster Variation 
Method (CVM) while efficient algorithms to include even infinite-range interactions 
in MC simulations have been devised |^ . MC simulations are also conceptually simple 
to understand, straightforward implement, and can be easily adapted to compute a wide 
variety of properties. 

Given these attractive features, it is not surprising that MC simulations of 
lattice models with cluster expansions determined from first-principles calculations 
have been employed in a variety of context for metallic, semi-conductor and ceramic 
systems, including the computation of: composition-temperature phase diagrams, 
thermodynamic properties of stable and metastable phases, short-range order in solid 
solutions, thermodynamic properties of planar defects (including surfaces or antiphase 
and interphase boundaries), and the morphology of precipitate microstructures [|6|, [f|, pi], 
221-0. 

With the steadily increasing availability of computational resources, MC 
simulations of lattice models are likely to become even more prevalent in computational 
material science modeling. Despite the technological advances in computational 
hardware, the human time needed to select the appropriate simulation parameters and 
to analyze the data has remained essentially constant over the years. We are now in 
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a situation where the main hmitation to the inclusion of atomistic MC simulations in 
standard thermodynamic software toolkits is not the lack of computational resources, 
but rather the difficulty of controlling a MC simulation in order to obtain the desired 
quantities. 

To address this problem, we have devised a variety of high-level algorithms that 
serve as an interface between the user and a traditional MC code. The user specifies 
the goals sought in a high-level form that our algorithm converts into elementary tasks 
to be performed by a standard MC code. For instance, our algorithms permit the 
determination of the free energy of a phase over all of its region of stability within a 
specified accuracy, without requiring any user intervention during the calculations. The 



proposed algorithms have been implemented in an easy-to-use software package [2(: 



In an effort to make this article self-contained, we first recall the basic theoretical 
results underlying semi-grand-canonical MC simulations, which are especially suited for 
determination of thermodynamic properties of alloys. We then proceed to describe the 
algorithms that enable the automated calculation of thermodynamic quantities from 
MC simulations. These algorithms address the following three questions, (i) How to 
determine when a MC simulation has run for a sufficiently long of time to provide the 
desired quantities within a given target precision? (ii) How to detect phase transitions? 
(iii) How to efficiently determine the composition-temperature phase boundaries of a 
phase diagram? We then provide a series of examples of applications that illustrate the 
features of our algorithms. 



2. Semi-Grand-Canonical Monte Carlo 



A very convenient way to obtain thermodynamic information regarding an alloy system 
is to perform MC simulations which sample a semi-grand-canonical ensemble, also known 
as the transmutation ensemble. In this ensemble, the energy and concentration of an 
alloy with a fixed total number of atoms (A^) are allowed to fiuctuate while temperature 
and chemical potentials are externally imposed. Although the results of this paper 
can be extended to general multicomponent systems, we are considering the case of 
binary alloys, for the benefit of notational simplicity. In an A — i? alloy, without loss of 
generality, only the composition x of element A and the difference in chemical potential 
between the two species n = ha — IJ^b needs to be specified. For conciseness, we simply 
refer to the quantity ^ as the "chemical potential" and x as the concentration. 

The natural thermodynamic potential (expressed per atom) associated with the 
semi-grand-canonical ensemble can be defined in terms of the partition function of the 
system as follows: 

/x) = In exp {-(3N {E, - /ix,))^ (1) 

where Ei and Xj are, respectively, the internal energy (per atom) and the concentration of 
state i, while (3 = 1/ {ksT) is the reciprocal of the temperature T and fc^ is Boltzmann's 
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constant. The thermodynamic potential is related to the Helmholtz free energy F 
through (p = F — fix. 

The potential can also be defined by the following total differential: 

d {(3(f)) = {E - fix) d(3 - (3x dfi. (2) 

where E and x are the system's average internal energy (per atom) and concentration. 
Equation (|^) enables the calculation of the thermodynamic function (p {j3, fi) through 
thermodynamic integration: 

= 0(/?o,/io) + "5 / {E - fix, -px) ■ d{p,fx), (3) 

where the integral is performed along a continuous path joining points (/?o, A*o) 
and that does not encounter a phase transition. Semi-grand-canonical MC 

simulations can be directly used to determine the values of E and x for any given (3 and 
fi. 

The potential at the initial point 0(/3o,/io) is usually chosen to be a convenient 
point where can be computed analytically. One possible starting point is the high 
temperature limit ||11|| , 

{(3, /i) - In (2) as (3^0 (4) 

where the limit is taken while keeping fi finite (implying that concentration x converges 
to 1/2 as /3 ^ 0). Another possible starting point is in the limit of low temperature at 
a chemical potential stabilizing a given ground state g. In this limit, a "single-spin flip" 
low temperature expansion (LTE) fll],^,^ can be used to obtain: 

(/3, fi) = Eg - fiXg - exp {-(3 {Aes,g - /iAr/,,^)) , (5) 

where Eg is the energy (per atom) of ground state g with composition x, while Asg^g is 
the variation in the system's total energy {NE) associated with changing the identity of 
the atom sitting at site s in ground state g, and Arjs^g is the variation in (Nx) associated 
with the same change. This infinite sum can be reduced to a finite number of terms by 
making use of the translational periodicity of the ground state. 

In order to avoid lengthy MC simulations, it is useful to possess a criterion to find 
the highest possible temperature where the LTE is accurate. One attractive possibility 
is to compute both from the LTE and from a mean field approximation (MFA). As 
long as both approaches give the same result within a user specified tolerance, we can 
be confident that "multiple spin fiip events" are rare and that the LTE has the desired 
level of precision. 

Using the LTE for each ground state of the system as a starting point to 
evaluate Equation (H), one can map out the whole potential surface 0" (/?, /i) associated 
with any given ordered phase a. Similarly, the high temperature limit can be 
used as a starting point to obtain 0" (/3, fi) for the disordered phase. This process 
is illustrated in Figure |l|. Once 0" {(3, fi) has been determined for all phases, the 
boundary between two given phases a and 7 can be located by identifying the locus 
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jya^y _ |(^^^) ; (J)"" (^jS , fi) = (p {(3,^)} where the potential surfaces intersect. At such 
intersections (/?, /i) G A"''*', the concentrations and x"^ of the two phases in equihbrium 
are simply given by the slopes of the potential surface at the point of intersection: 



The calculations of thermodynamic quantities based on semi-grand-canonical MC 
simulations offer several advantages over their canonical counterparts. 

(i) For a given value of the control variables /i), the thermodynamic equilibrium 
of the system is never a phase-separated mixture, implying that the calculated 
quantities always reflect the properties of pure phases, free of interfacial 
contributions that would otherwise not be negligible due to the finite size of the 
simulation cell. (Note that the MC simulation cell must be commensurate with the 
periodicity of the equilibrium structure for this property to hold.) 

(ii) The potential can be very directly determined through a simple thermodynamic 
integration procedure (Equation (^)) where each required quantity takes the form 
of a thermodynamic average of computationally inexpensive quantities {E, which 
must already be computed in order to perform the simulation, and x, which is 
trivial to obtain). 

(iii) Phase boundaries can be located by looking for intersections between curves, 
a criterion which is somewhat simpler to implement than the common tangent 
construction (although both methods are formally equivalent). 



Self-driven Monte Carlo simulations 



6 



3. Algorithms 

While the methods described in the previous section in principle enable the 
determination of thermodynamic quantities from MC simulations, a number of practical 
issues need to be addressed before this process can be fully automated. 

A typical MC code requires a substantial amount of user intervention and relies 
on the user's physical intuition to identify when a simulation is successful and reliable. 
The type of user interventions needed include: 

(i) Providing a variety of input parameters controlling the computational accuracy 
that need to be determined by trial and error. 

(ii) Requiring the user to monitor the simulation periodically to ensure that the system 
has not undergone an unexpected phase transformation. 

(iii) Performing a substantial amount of post-processing in order to obtain a phase 
diagram. 

In this section we introduce algorithms to perform these tasks automatically and 
provide a formal justification of their applicability. 

3.1. Reaching a target precision 

Three parameters control the precision of a thermodynamic quantity computed through 
MC simulations: the simulation cell size, the equilibration time and the averaging time. 
While the analysis of the effect of simulation cell size on the precision is beyond the 
scope of the present paper, this section provides a sound basis for the selection of the 
two remaining parameters. The equilibration time is the time the system takes to 
reach thermodynamic equilibrium from a given initial configuration. Once equilibrium 
is reached, the instantaneous value of a given quantity (e.g. the energy) then needs to 
be averaged over a certain number of MC steps in order for its average to be sufficiently 
accurate. This defines the averaging time, which is the first parameter we will focus on. 

3.1.1. Averaging time While the approach presented in this section relies on standard 
statistical results, it is nevertheless included for the purpose of introducing our notation 
and clarifying the assumptions made. Consider a microscopic quantity taking the value 
Qt at step t and whose expectation value is equal to the desired thermodynamic quantity 
Q. If L observations of Qt are available, a natural estimator of Q is the average: 



The precision of this statistical quantity can be quantified by calculating its variance. 



_ 1 ^ 



(8) 



In general, the variance of Q[i^l] given by: 




t=l s=l 



(9) 
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where V\t-s\ is the covariance between Qt and Qg- The fact that the sequence Qt is 
stationary when the system has reached thermodynamic equihbrium enables a simple 
estimator Vi of the covariance structure Vi of the sequence through the equation 

Vi = j^Y.QtQw-Ql.L]- (10) 
^ * t=i 

While Equation ^ is very general, it is also very computationally intensive, requiring 
on the order of operations to estimate the variance of Q^i Following most of the 
literature on MC simulation, we consider the more computationally tractable special 
case obtained when the covariance structure Vi is assumed to be of the form 

Vi = Vp-\'\. (11) 

This approximation can be motivated as follows. A Markov chain can be represented 
by a linear superposition of components with exponentially decaying autocorrelation 
function. Each component is characterized by a certain correlation length and the 
component with the longest correlation length will typically give by far the largest 
contribution to the variance. Considering only the most slowly decaying component 
should thus provide a good approximation to the true variance of the quantity of interest. 

Under the assumption that the averaging time is much larger than the longest 
correlation time,| "boundary" effects can be neglected and Equation then reduces 
to 

-| L oo 

{Qhm) = E (12) 

t=i i=— 00 

Evaluating the geometric sum then yields: 



Var (Q 




1-./ 

Note that, since p = exp (— 1/r) where r is the characteristic relaxation time. Equation 
(|l3|) can be written, in the limit of large r, in the more familiar form Var {Q[i l]) = 
2tV/L. 

While V can be directly estimated from Equation (|10|), setting / = 0, the parameter 
p can be found in a variety of ways. We found the following algorithm to be both fast and 
accurate: Find the time interval I* such that Vi*/Vo = | and then set p = 2~^/'*. This 
approach offers the important advantage that it uses, as an input, only the correlation 
between relatively distant values of Qt- This ensures that p will be determined by the 
most slowly decaying components of the Markov process, which are the ones that the 
functional form l^p''' was meant to model. It is also a very fast algorithm, since a simple 
bracketing algorithm allows the determination of p in about L log L operations. In the 
event that there are more than one /* satisfying the equality, we chose the smallest one. 
The rationale behind this choice is that there is a small but nonzero chance that the 

I This assumption should not be a concern because (i) this assumption tends to be violated only for 
averaging times such that the desired target precision on Q has not yet been reached, and (ii) our 
expression overestimates the true variance when this assumption is violated. 
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estimated covariances are such that Vi/Vq = 1/2 although the true covariances are not 
such that Vi/V = 1/2. Choosing the smallest I* makes the procedure less sensitive to 
these spurious events since the probability of occurence of one such event increases with 
the range of values of / considered. 

In short, the averaging time L needed to reach a given precision can be found by 
computing the variance of Q[i^l] through Equation (|T3|) periodically and stopping the 
simulation when the variance goes below a user specified threshold: 

^ar(Q[i,^])<(^)', (14) 

where p is the target precision on the value of Q while 

2„ = \/2erf"^(l -a) (15) 

is the critical value associated with the desired level of confidence a. For instance, if the 
error on Q is required to be less than 10~^ with a probability of 99%, then p = 10~^, 
a = 0.01 and 2„ = 2.576. 



3.1.2. Equilibration time When devising a criterion to determine whether a simulation 
has reached thermodynamic equilibrium, approaches based on identifying correlation 
lengths can be unreliable because the system may be very far from equilibrium 
during equilibration. For this reason, we employ a relatively simple, and yet very 
robust approach that is based on the main defining characteristic of thermodynamic 
equilibrium. If the average of the quantity Qt between steps ti and t2 is not statistically 
significantly different from the average of Qt between steps ^2 and ta, then the hypothesis 
that equilibrium has been reached at step ti and beyond cannot be rejected. Of course, 
the power of this test increases as ^2 — ti and t^ — t2 increase and this definition of 
equilibrium thus depends on the target precision p we seek to reach. This is to be 
expected since, strictly speaking, thermodynamic equilibrium is never reached in a 
simulation of a finite length. Any definition of equilibrium must include an arbitrary 
cutoff and it is natural to choose this cutoff to be a function of the precision we require. 

Both the equilibration and averaging time needed to reach a precision p can thus 
be determined through an algorithm that can be outlined as follows (refer to Figure 0). 
Let L denote the number of samples of Qt collected so far. Consider a range of values 
of ti G [1, L] and for each, check whether 

AQ = Q[t^^(ti+L)/2] ~ Q[{ti+L)/2,L] < (16) 

where Q[ti,t2] denotes the average of Qt for t G [^1,^2]- Then, verify that Var (Qj^^ j^]) 
satisfies Equation ([Mf) . If both conditions are satisfied, then Q^^^ ^] provides an estimate 
of Q with a precision p and with a confidence level a. 

In its simplest form, the implementation of this algorithm involves storing all 
samples of Qt collected so far and requires a loop over values of ti which make the 
complexity of the algorithm 0{L^). Various techniques help reduce the resources 
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a) 





Figure 2. Criterion for testing whether thermodynamic equiHbrium has been reached 
in a MC simulation. In part a), the two blocks of data have significantly different 
means, implying that the hypothesis that equilibrium has been reached at ti is rejected. 
In part b), the two blocks of data do have similar means, indicating that the hypothesis 
that equilibrium has been reached at t'l cannot be rejected. 



required. First, sampling Qt every MC pass,|| instead of at every spin flip, results 
in essentially no loss in precision since the correlation length of the simulation typically 
extends much beyond a MC pass. Second, one can limit the search for ti to values 
that are integer multiple of some block size that is allowed to increase linearly with 
the total number of MC passes L. In this fashion, the number of operations needed to 
determine ti does not increase as the simulation progresses. The value of ti thus found 
by a search constrained to block boundaries differs at most by Lh from the ti which 
would have been found without block constraints. In the worst case, the simulation 
may have to be run for an additional Lf, passes and it follows that the length of the 
simulation will be at most (1 + Lb/L) larger than the optimal computational time that 
would have been needed without block constraints. A computationally convenient way 
to implement this scheme is to multiply the block size by two whenever the total 
number of passes L doubles. Thanks to this approach, the entire history of values of Qt 
does not need to be stored, only block averages of Q, V and p are needed, as illustrated 
in Figure |^. When the block size doubles, the values of Q, V and p for two consecutive 
blocks can be combined to obtain the corresponding value for the new block having 
twice the length.[j]] 

§ A MC pass consists of a number of attempted spin flips equal to the number of lattice sites in the 
simulation cell. 

|j The only drawback is that the resulting estimate of p does not make use of all the available data 




T'i, = 2Tt 

Figure 3. Block-wise storage scheme, a) Only the mean variance and correlations 

within blocks of width Lf, needs to be stored (the mean is represented by a thick 
horizontal line), b) As the simulation runs, blocks are periodically merged so that 
storage requirements do not grow. 



3.2. Detecting phase transitions 

Phase transitions can be identified by locating singularities in some function Q{C), 
where Q is some thermodynamic quantity (such as energy, concentration or an order 
parameter) and C is some control variable (such as temperature or chemical potential) . 
In this section, we propose a simple algorithm to perform the detection of such 
singularities while accounting for the fact that quantities obtained from MC are 
intrinsically noisy. Our algorithm is applicable to the detection of both first-order and 
second-order transitions, although the power of the test is clearly higher for first-order 
transitions. 

An analytic function has the property that the knowledge of its shape in some finite 
interval enables the prediction of its value outside of that interval. The failure of this 
extrapolation procedure in the presence of a singularity permits its detection. Formally, 
an analytic function can be extrapolated using a Taylor series. We approximate this 
operation by fitting the output of the MC simulation by a polynomial. 

Consider n + 1 estimates Q^, . . . , Qn+i of the quantity Q (C) obtained for a sequence 
of values of the control variable Ci, . . . ,Cn+i- The general principle behind our 
singularity detection procedure is to fit a polynomial through the first n points and 

points, because the product QtQs for t and s belonging to different blocks cannot be determined from 
within-block averages. Although it is inefficient, our estimate of p is nevertheless unbiased and the loss 
of efficiency is negligible if the block size is large relative to the correlation time. 
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Figure 4. Criterion for detecting a phase transition. If the sequence (Q„, C„) predicts 
a value of (^Q„^i, C'n+ij but the actual value, {Qn+i, C„+i) is statistically significantly 
different, then a phase transition must have occurred. Both the noise in the prediction 
(represented by dashed lines) and the noise in Q^+i (represented by an error bar) must 
be taken into account to determine statistical significance. 




if the observed value Q^+i is far from the value predicted from the polynomial fitted to 
the previous data, a singularity must be present between point n and n + 1 (see Figure 
To make this algorithm practical, two questions need to be answered. First, how 
do we choose the order of the polynomial to be fitted? Second, how do we select the 
maximum allowable prediction error that can be tolerated without claiming the presence 
of a singularity? 

To answer the first question, we make use of the well-known cross-validation 
criterion |^,^: Choose the order k of the polynomial that minimizes the quantity: 

cv=-j:{Q,-Q:y (17) 

^ 1=1 

where Q* is the value of predicted from a polynomial least-squares fit to the remaining 
n — 1 points: Qi, • • • , Qi-i, Qi+i, ■ ■ ■ , Qn- Intuitively, this criterion tests the predictive 
power of the polynomial using the points Qi, . . . before we attempt to use it for 
the prediction of Qn+i- Once the order k of the polynomial has been chosen, the 
standard theory of least-squares provides the distribution of the predicted value 
Qn+ii which can be used to determine whether the prediction error is statistically 
significantly different from zero. For i = 1, . . . , n and j = 1, . . . , A; + 1, let 

x.,= {c,y-' (18) 

Y^ =Q„ (19) 

so that the equation defining the least-squares problem can be written in matrix form 
as y = Xa + e, where e is an n x 1 vector of error terms. The vector a of the coefficients 
of the polynomial can be estimated by 

^1 



X'Y (20) 
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while the covariance matrix of a is given by 

V = a'{x^Xy' (21) 

where is the variance of the residuals Ei of the fit. The quantity cr^ can be estimated in 
two asymptotically equivalent ways. Equation (0) provides an estimate of the variance 
of any of the Q^, which can be averaged over i to yield a single, more accurate estimate: 

a' = -f: Var (Q,) . (22) 

^ i=l 

Alternatively, the residuals of the least-squares fit can be used: 

= \\Y -Xaf / (n-k-l) . (23) 

While Equation (^2]) is more accurate (especially if n is small), because it makes use 
of additional information. Equation ( p3D is simpler to implement as it minimizes the 
amount of information that needs to be transferred between the Monte Carlo simulation 
code and the phase transition detection routine. 

The predicted value of Q for point + 1 is then given by 

g:+i = x^a (24) 

where Xj = [Cn+iY''^ for j = 1, . . . , + 1, while the variance of this prediction is given 
by: 

V = x^Vx. (25) 

If there is no singularity, the variance of the difference between Q^+i and Q*n+i is 
v + a"^ because (i) the variance of Q^+i is v while the variance of Qn+i is ^"^i (ii) Qn+i 
Qnj^i are independent. An asymptotically valid statistical test thus consists in rejecting 
the hypothesis that there is no phase transition if the prediction error Qn+i — Qj^+i 
such that^ 

Qn+l-Qn+l\>Zay/VT^. (26) 

For Za defined as in Equation ([T5|) , the probability that a phase transition is incorrectly 
identified is less than a. The relevant quantities are represented in Figure ^. 

The Equations presented above rely on the assumption that the deviations away 
from the true thermodynamic quantity, (^Q^ — Q {Cij^ for z = 1, . . . , n + 1, are 
statistically independent. This is an appropriate assumption if the system is given 
sufficient time to equilibrate every time the value of the control variable C changes. 
Our method also assumes that the deviations (^Q.i — Q (Cj)^ all have the same variance. 
Both of these requirements are automatically satisfied if the points Qi, . • • ,Qn+i ^-re 



is 



generated using the algorithm presented in Section 3.1 



^ Note that Qn+i must be an unbiased predictor of Q (C„+i) for this test to be vaUd, which is not the 
case for finite n due to the finite order of the approximating polynomial. Fortunately, the bias (squared) 
goes to zero as n — > oo and eventually becomes negligible relative to tr^ , ensuring the asymptotic validity 
of the procedure. Note that the variance v of the prediction is also negligible asymptotically relative 
to cr^. We nevertheless include the variance term v in Equation (|2^) in order to improve the accuracy 
of the test for finite n, since we do have an estimator of the variance (unlike the bias). 
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For the sole purpose of locating phase transitions, one particular type of 
thermodynamic function enables an especially sensitive test for phase transitions: order 
parameters. Unfortunately, order parameters that are invariant under any possible 
symmetry operation of the lattice are often expensive to compute. However, if one 
always traverses phase transitions from an ordered phase toward a disordered phase (or 
another ordered phase), the task is dramatically simplified. When a MC cell is initialized 
with one specific ground state (i.e. an ordered phase) and before a phase transition is 
reached, there is a vanishing probability that the system will fluctuate enough to reach a 
state that is different but symmetrically equivalent to the original state. Thus, an order 
parameter can be computed by simply keeping track of the number of sites that host an 
atomic specie that is different from the specie found at the same site at initialization. 
There is no need to check for every possible symmetry-related ordered configuration. 

The algorithms just presented are especially useful when one desires to construct 
the potential surface (/5, fi) of a phase (using thermodynamic integration over the 
paths depicted in Figure |I]) without requiring the user to manually determine when to 
stop the thermodynamic integration along a given direction because a phase transition 
has occurred. The automated construction of such free energy surfaces provides a very 
natural pathway through which an atomistic MC simulation can be used to provide 
input to CALPHAD-type calculations based on phenomenological free energy models. 

3.3. Phase boundary tracing 

When one is solely interested in determining a system's phase diagram, it would be 
computationally advantageous to be able to follow the boundary of a phase without 
first calculating the potential over the whole region where this phase is stable. In this 
section, we describe how, in the case of a first-order transition, the entire boundary of a 
two-phase equilibrium can be determined, starting from a single point where a transition 
is known to occur. We first focus on a simple method that neglects the unavoidable 
presence of statistical noise in MC simulations before describing how the presence of 
this noise can be handled. 

Consider a value of the reciprocal temperature (3 and of the chemical potential fj, 
that is known to stabilize a phase-separated mixture between phases a and 7. The 
knowledge of such a pair (/3, /i) could be provided, for instance, by the phase transition 
detection algorithm presented in the previous section, assuming that the hysteresis 
loop of the phase transition is sufficiently narrow (which tends to be case at high 
temperature). Alternatively, if T = 0, the chemical potential that simultaneously 
stabilizes two ground states a and 7 can be found by a simple common tangent 
construction from the knowledge of the energy and composition of each ground state. 
At such a two-phase equilibrium point, the potentials of each phase are equal. Our 
task is to find, as reciprocal temperature /3 changes by an amount d(3, the changes in 
chemical potential needed to preserve the equality between each phase's thermodynamic 
potential. By iterating this process, one can define a whole path /i (/5) that traces out 
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the value of the chemical potential at the phase transition as temperature changes, and 
thus find the concentrations Xa and of each phase in equilibrium as a function of 
temperature. The differential equation whose solution provides the desired path fi 



can be found by taking the total differential of the equation 
the identity between the potentials (p of each phase: 



P following from 



dp 



dp 



d(3 + 



dfi 



{E" - /ix") d(3 - 
dix (^^ - E") /i 
dp 



{E'< - /ix^) dp - px^dfi 



(27) 

(28) 
(29) 



/3(x^-a;") p' 

where E"' and are the energy and concentration of phase a at reciprocal temperature 
P and chemical potential fi (and similarly for phase 7). All the quantities needed can 
be directly obtained from MC simulations. Thanks to the fact that both phases are 
metastable for values of the chemical potential in the neighborhood of the true phase 
transition, it is possible to obtain E and x for both phases at the same chemical potential. 
A very efficient way to implement this algorithm is to keep two simulations in memory 
at the same time — one for the phase a and one for the phase 7. In this fashion, the 
final configuration of the simulation of phase a at reciprocal temperature P can be used 
as an initial configuration for the simulation of phase a at temperature P + dp (and 
similarly for phase 7). The long equilibration time associated with the crossing of a 
phase transition is thus avoided. Interestingly, this algorithm can be directly used even 
when the two phases in equilibrium are based on a completely different parent lattice. 

We now address the issue that the quantities obtained from MC are contaminated 
by noise. The presence of noise causes the solution to Equation (p9D calculated from 
MC to deviate from the true path fi{P). This eventuality brings up two questions. 
First, is the solution to Equation ( pOj) stable or unstable with respect to the presence 
of small perturbations? Second, what can one do if the noise becomes so large that the 
calculated chemical potential leaves the region where both phases are metastable? 

To answer the first question, consider a perturbation e (P) of the chemical potential 
away from its value /i (/?) at the phase transition at reciprocal temperature p. We now 
calculate how this error is further propagated when solving Equation To the first 
order, the difference between the potential of the two phases is P (x'^ (P) — x" (P)) e (P). 
As temperature changes. Equation (|29|) updates fi so as to preserve this difference. We 
thus have the equality: 

(x^ iP) - x" (/?)) e iP) = (x^ iP + dp) -x"{P + dp)) e{P + dp) . (30) 

Taking the limit of infinitesimal dp and rearranging yields: 
de d 



(31) 



and it follows that the small perturbation e is unstable as P increases when the 



difference in concentration between the two phases in equilibrium 



X ' 



X 



decreases 



with increasing p. The practical consequence of this observation is that integrating 
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Figure 5. Criterion for detecting when the phase boundary tracing algorithm has 
exited the region of simultaneous metastability of the two phases in presence. The 
point [x"^i, Pn+i) is better predicted by {x^*^i,Pn) than by [x"*^i, Pn) ■ Hence, the 
MC simulation of phase a must have undergone a transition to phase 7. 



Equation (^91) from high temperature to low temperature yields a method more robust 
to statistical noise, because \x'^ — x°'\ usually increases with decreasing temperature T. 

We now address the second robustness issue, namely, what can we do if the 
algorithm wanders outside of the region of metastability of either phase? It is 
straightforward to detect when this happens by simply checking whether {x"^ — x°'\ is 
not statistically significantly different from zero. One then needs to determine which 
of the two MC simulations (for phase a or 7) underwent a phase transition. This can 
be accomplished by keeping track of the concentration of each phase (x", . . . , and 
x1, . . . ,xX) at temperatures previously visited . . . and using them the verify 
whether the new concentration x^_^_i ~ at Pn+i is best predicted by the sequence 
or x] , using the extrapolation described in the previous section (see Figure]^). The 
sequence that gives the worst prediction will point to the MC simulation that underwent 
a phase transition. 

Without loss of generality consider the case when x" < x'^ and when the simulation 
of phase a transformed to phase 7. We can then "recenter" the chemical potential to a 
good estimate of its true value at the phase transition as follows (refer to Figure |^). 

(i) Gradually decrease the chemical potential of the MC simulation of phase a until 
its concentration can again be predicted by the sequence x" within the magnitude 
of the statistical noise. Let /i denote the value of the chemical potential when this 
happens and save the state S of simulation cell at that point. 

(ii) Then, gradually increase the chemical potential of the MC simulation of phase a 
until its concentration can be predicted by the sequence x] within the magnitude 
of the statistical noise. Let JI denote the value of the chemical potential when this 



Self-driven Monte Carlo simulations 



16 




Figure 6. Algorithm to locate the center of the region of simultaneous metastability 
of the two phases in presence. In step 1, the chemical potential in the simulation 
of phase a is gradually decreased until concentration is once again best predicted by 
xf* . In step 2, the chemical potential is increased until the simulation transform back 
to phase 7. The two extreme values of the chemical potential /i and JI thus defined 
bracket the true value of /i at the phase transition. 



happens. 

(iii) Restore the state of the simulation cell of phase a to the saved state S and set the 
chemical potential to /i = (jl + jjj /2. 

(iv) Rerun the simulation for both phases a and 7 at chemical potential /i and continue 
the integration Equation (|29|) with the values of E and x obtained for each phase. 

This approach results in a very robust algorithm that can run without user 
intervention to map out the whole two-phase equilibrium. Thanks to this safety 
mechanism, we found that it is feasible to integrate Equation (^9[) towards increasing 
temperature, despite the potential for instability with respect to noise. This proves 
useful since it is very easy to find a starting point for our algorithm at T = 0. The 
boundary tracing algorithm may have to recenter itself a few times on its way to the 
end of the two-phase equilibrium. But once the location of the phase boundary at high 
temperature is known, a smooth phase boundary can then be calculated, starting from 
the highest temperature and following the phase boundary in the direction of decreasing 
temperature. 

The end of the two-phase region can be detected by monitoring the occurrence of 
either one the following events (see Figure 0): 

(i) The concentrations and x'^ are statistically significantly different, but either 
or x^ is not predicted by the corresponding sequence xf or xjwithin statistical 
uncertainty. 
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Figure 7. Criterion for detecting tlie end of a two-pliase equilibrium. In Case 1, a 
new phase S has appeared, as indicated by the fact that x'^_^,i lies outside of the error 
bar of the prediction of concentration for either phase a or 7. In Case 2, a congruent 
or critical point has been reached since the concentrations of the two phases are no 
longer statistically significantly different. 



(ii) The concentrations and x'^ are not statistically significantly different, but 
the sequences xf and x] both predict the same concentration, within statistical 
uncertainty. 

Case 1 indicates that a third phase S has become stable (i.e. a eutectic or peritectic 
point), while Case 2 indicates that the phase boundaries have met at a common 
composition (e.g. at a congruent point or a critical point). In the first case, the 
algorithm will actually continue slightly beyond the true eutectic or peritectic point 
because of the presence of hysteresis at a first-order phase transition. However, the true 
location of the three-phase equilibrium can be determined from the intersection of the 
a — 6 and 6 — 'j phase boundaries, once they have been determined. 



4. Applications 

This section presents a number of examples of automated calculations of thermodynamic 
properties using the previously introduced algorithms. In each example, the cluster 
expansion describing the energetics of the alloy system was obtained from a fit to first- 
principles structural energy calculations using the MAPS ||3^ code, which automates 



the process. A description of the algorithms enabling the automatic construction of 
such a cluster expansion, as well as a description of the interaction parameters defining 
each cluster expansion used below, can be found in [p3| . 

Our first example is the determination of the complete potential surface (T, fi) 
of the solid-state phases of the MgO-CaO pseudo-binary alloy system (see Figure 
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iiJ. SKIMI 

Figure 8. Potential surface for the MgO-CaO pseudobinary alloy system. 



g). As the MgO-CaO system is phase separating (as shown in Figure |T2D, the 
thermodynamic potential surface exhibits a characteristic butterfly shape. Below the 
critical temperature, the surface intersects with itself at a nonzero angle, indicating the 
presence of a miscibility gap. Since the MC simulation is able to sample metastable 
phases, for some values of the chemical potential, one can determine (T, /i) for 
both phases in the region of metastability, which explains the double-valued nature 
of the potential surface at low temperatures. Above the critical temperature, (T, /i) 
becomes single- valued and varies smoothly as fi varies. (Note that, in reality, the 
alloy melts before reaching that point at atmospheric pressure.) Figure |^ shows 
the GibbsQ free energy surface G {x,T) = (p {T, fi {T,x)) + fi{T,x)x for the same 
system and, in this more familiar representation, the presence of miscibility gap at 
low temperature is clearly recognizable. The calculated Gibbs free energy surface 
can be directly used to fit standard solution models, such as an ideal solution model, 
plus a polynomial in temperature and concentration to account for non-ideality. The 
order of the polynomial needed to represent the MC output can, for instance, be 
determined through cross-validation, as defined in Section p.2| . The resulting very 
compact representation of the thermodynamics of the alloy system can be used in 
CALPHAD-type calculations to supplement existing experimental data or to fill in data 
not yet determined experimentally. 

Figure 10 shows the result of a similar calculation in the case of a ordering alloy 
system, the Ti-Al system. For clarity, the metastable portions of the (T, /i) surface are 
omitted so that the potential is single-valued everywhere. The order-disorder transition 
is clearly visible from the kink in the surface, which is highlighted by a thick curve. The 

+ The difference between the Gibbs and Helmholtz free energies can be neglected at atmospheric 
pressure in solid state systems because the product of pressure P with changes in specific volume {AV) 
are small. 
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Figure 9. Gibbs free energy for the MgO-CaO pseudobinary alloy system. Common 
tangent constructions are shown at selected temperatures. 




corresponding Gibbs free energy representation of the same data is shown in Figure |Tl] . 

Figure 12 illustrates the usefulness of the boundary tracing algorithm described 
in Section |3.3| . The examples shown even include a case (Ti-Al) where a equilibrium 
between two phases differing by their lattice type: the LIq is fcc-based while the DOig 
and the Ti solid solution are hep-based). These examples show that one can obtain phase 
diagrams from first-principles without first obtaining the complete (T, /z) surface. 

It is important to note that, in these examples, the MC code requires very little 
intervention or input from the user. Thanks to the algorithms introduced in Section 

the code is able to autonomously determine, at each value of the control variables 
(T, /i) , when equilibrium has been reached and how long the simulation needs to be run 
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Figure 12. Sample phase diagrams calculated with the automatic boundary tracing 
algorithm. 
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in order to reach the target precision. The code is also able to determine when the 
region of stability of one phase has been exited, thanks to the algorithms described in 
Section p?2| , so that the user does not need to specify when to stop the thermodynamic 
integration. The only input parameters the user needs to specify are 

(i) which phase to scan (the code can determine the value of the chemical potential n 
stabilizing the desired two-phase equilibrium at OK); 

(ii) the target precision (here, 0.1 at % precision on the concentration variable is used 
as a criterion); 

(iii) the size of the simulation cell (in the present study, the simulation cell is always 
chosen such that it contains a sphere of a diameter equal to at least 10 times the 
nearest-neighbor interatomic spacing). 

5. Conclusion 

The algorithms we have introduced enable researchers to focus on high-level concepts 
when employing MC simulations. All the tasks that make traditional Monte Carlo 
simulations so tedious have been formalized into algorithms that can autonomously 
initialize and control the appropriate MC simulation. The user only needs to specify 
high-level goals, such as the target precision or which phase to focus on. This automation 
of MC simulations, along with the earlier work that enabled the automation of cluster 
expansion construction signifies that it is now possible for large community of 
researchers to employ first-principles calculations to augment and complete existing 
experimental thermodynamic databases without facing the steep learning curve that 
has traditionally been associated with first-principles thermodynamical calculations. 
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